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Abstract: We present the cosmic shear signal predicted by two viable cosmological models 
in the framework of modified-action f(R) theories. We use f(R) models where the current 
accelerated expansion of the Universe is a direct consequence of the modified gravitational 
Lagrangian rather than Dark Energy (DE), either in the form of vacuum energy /cosmological 
constant or of a dynamical scalar field (e.g. quintessence). We choose Starobinsky's (St) and 
Hu & Sawicki's (HS) f(R) models, which are carefully designed to pass the Solar System 
gravity tests. In order to further support - or rule out - f(R) theories as alternative candidates 
to the DE hypothesis, we exploit the power of weak gravitational lensing, specifically of cosmic 
shear. We calculate the tomographic shear matrix as it would be measured by the upcoming 
ESA Cosmic Vision Euclid satellite. We find that in the St model the cosmic shear signal is 
almost completely degenerate with ACDM, but it is easily distinguishable in the HS model. 
Moreover, we compute the corresponding Fisher matrix for both the St and HS models, thus 
obtaining forecasts for their cosmological parameters. Finally, we show that the Bayes factor 
for cosmic shear will definitely favour the HS model over ACDM if Euclid measures a value 
larger than ~ 0.02 for the extra HS parameter tt-hs- 
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1. Introduction 

It is now widely accepted that a number of different cosmological observations points out a 
lack in our current understanding of the Universe. In particular, the temperature anisotropy 
spectrum of the Cosmic Microwave Background (CMB) [1, 2, 3, 4, 5, 6, 7], the Hubble 
diagram of Type la Supernova? (SNela) [8, 9, 10, 11, 12, 13, 14, 15, 16] and the clustering 
properties probing the cosmic Large-Scale Structure (LSS) [17, 18] are concordant pieces 
of evidence in favour of a present period of accelerated expansion of the Universe. In the 
current cosmological model A Cold Dark Matter (ACDM), this observational evidence and 
the problem of the missing mass in the dynamics of galaxies and galaxy clusters and in the 
LSS of the Universe [19, 20, 21, 22, 23, 24] are excellently reproduced and solved by assuming 
a spatially flat Universe dominated by cold Dark Matter (DM) and a vacuum Dark Energy 
(DE) in form of a cosmological constant A. 

However, this scenario has serious drawbacks: if A is interpreted as vacuum energy, its 
value is 120 orders of magnitude smaller than what is expected from quantum field theory; 
the coincidence and fine-tuning problems do not seem to have any natural explanation and 
DE has to represent ~ 70% of the total energy budget of the Universe. This circumstance has 
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motivated the search for alternative solutions. They mostly rely on scalar fields with suitable 
potentials which provide a varying A term or other DE components or even unified DM and 
DE fluids [25, 26, 27, 28, 29, 30, 31] with exotic properties. 

Recently, a different route has started to be followed. The cosmological constant was 
originally proposed as a geometric term, i.e. a constant term in the left-hand (geometric) side 
of the Einstein equation, thus adding no extra component in the stress-energy tensor of the 
fluid filling the Universe. By generalising this approach, one can argue whether it is possible to 
reproduce the current cosmic accelerated expansion by adding a non-constant time-dependent 
term in Einstein's tensor. The effort of modifying and generalising the Hilbert-Einstein action 
of gravity actually dates back to just few years after Einstein's seminal papers (e.g. [32] for 
a historical review), and it has been also proposed by Starobinsky [33] in order to explain 
the cosmic inflation in the early Universe. This idea has been suggested again nowadays for 
the purpose of correctly describing the current accelerated expansion of the Universe without 
any exotic fluid [34, 35, 36, 37, 38, 39]. These modified-action theories of gravity are widely 
known as f(R) theories, because in the gravity Lagrangian the Ricci scalar R is replaced by 
a generic function f(R). 

However, GR is a well-tested theoretical framework, at least with respect to the Solar 
System scale of distances. Therefore, any f(R) theory which attempts to solve the late-time 
acceleration problem has to face the Solar System tests of gravity. Recently, two models 
carefully designed to pass the local gravity tests but still providing an accelerated cosmic 
expansion have been proposed [40, 41]. 

In this paper we use weak gravitational lensing to test and constrain these models. Indeed, 
weak lensing is a powerful tool [42, 43, 44, 45, 46, 47, 48, 49, 50]: gravitational lens effects are 
due to the deflection of light occuring when photons travel near matter, i.e. in the presence of a 
non-negligible gravitational field. The cosmic convergence and shear encapsulate information 
about both the source emitting light and the structures that photons cross before arriving at 
the telescope. Therefore, weak lensing allows to explore both the basis of the cosmological 
model and LSS of the Universe, in other words it brings information about the geometry and 
the dynamics. Hence, the study of the power spectrum of weak lensing can be a crucial test. 

The structure of this work is as follows. In § 2, we present the theory of f(R) models, both 
at background level and in the linear theory of cosmological perturbations. In § 3, we recast 
the most important observables of weak gravitational lensing, the convergence and the cosmic 
shear, and their power spectra. Then, we also outline how to construct the tomographic shear 
matrix. § 4 illustrates the cosmic shear Fisher matrix, whilst §4.1 the Bayes factor and its use 
in the framework of model selection. Finally, we show our results in § 5. In § 6, conclusions 
are drawn. 

2. f(R) Cosmology 

In modified-action theories of gravity, the standard Hilbert-Einstein Lagrangian density of 
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GR in the ACDM model 

is replaced by a more general function f(R) of the Ricci scalar R (see [51] for an exhaustive 
review) , 

Sf = I&. (2.2) 

Modified Einstein's equations are given by the variation of the Lagrangian with respect 
to the metric. This gives, after some manipulations and modulo surface terms, 1 

f, R R^ - V^f,R + (of tR - M 9ilv = SirGT^, (2.3) 

where t R denotes a derivative with respect to R, R^ v is the Ricci tensor, g^, the spacetime 
metric and T^ u the energy-momentum tensor. 

In these modified theories of gravity, the usual massless spin-2 graviton is not the only 
carrier of the gravitational interaction. Indeed, there is also a scalar degree of freedom - often 
dubbed scalaron -, which can be conveniently described by the function 4> = — 1. The 
evolution of the field cf> obeys to the trace of Eq. (2.3) 

where the potential V is related to R by 

^ = ±(2f-Rf tR )f, RR , (2.5) 

and T is the trace of the stress-energy tensor. 

Albeit the choice of the functional form of f(R) is in principle completely free, there 
are some constraints which help to avoid such an arbitrariness. First, modifying the gravity 
Lagrangian leads to deviations from GR at all scales. Particularly, the presence of the scalar 
degree of freedom 4>, which couples to matter, is responsible for a long-range fifth force. 
At the scale of the Solar System, this new force can possibly lead to wrong values of the 
PPN parameters [52, 53, 54]. A usual way to solve similar problems is to reassociate high 
densities with high curvatures, so that <p becomes very massive and the fifth force escapes 
any detection. Such a behaviour of the model, namely a large and positive mass squared 
term m 2 = (8itG) 2 p/3 in high curvature environments, is usually called "chameleon effect." 
On the other hand, the Big-Bang nucleosynthesis and early Universe data, for example the 
temperature anisotropies of the CMB, strictly constrain the choice of the cosmological model. 
In other words, one wants to recover GR at high redshift z. Similarly, the background 
evolution of the cosmos should not be too different from what predicted by the ACDM model, 
since this model reproduces SNela data excellently. 



x We use units such that c = 1 and signature { — , +, +, +}, where Greek indeces run over spacetime dimen- 
sion, whereas Latin indeces label spatial coordinates. 
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The constraints which must be satisfied by f(R) can be summarised as 



hm/(i?) = 0, (2.6) 
lim f(R) = const., (2.7) 

R— s-oo 

f,R\ R » m * > °- ( 2 - 8 ) 

A further condition which ensures that the solution is stable at high curvatures should be 
included. This translates into 

f,RR\ R>>m 2 > 0. (2.9) 

Among the possible choices left out by the above conditions, in this paper we will consider 
two popular classes of f(R) models which we briefly describe in the following, namely the 
models proposed by [40] and [41]. 

Starobinsky Starobinsky's model [40] (St, hereafter) is described by the function 



f(R) = R + XR, 



2 \ ~ n St 



R 2 , 

1 + -J 



R* 



(2.10) 



with R+ a scaling curvature parameter and A and nst two positive dimensionless constants. 
It is worth noting that there is no cosmological constant, because /(0) = 0. Nonetheless, 
in high curvature regime, one recovers an effective A-like term, since f{R S> R±) — R — 
2A( cff ), where A( eff ) = XR+/2. Moreover, in the early Universe the Hubble parameter H is 
the same as in a ACDM model with an effective matter fraction fi^f^ = 1 — \R±/ {QHq 2 ), 
where Hq = 100/ikm s _1 Mpc -1 is the Hubble constant. This guarantees that the Big-Bang 
nucleosynthesis constraints are satisfied. 

Hu & Sawicki Hu k, Sawicki [41] proposed an f{R) model (HS, hereafter) described by 

f w = R - m i + c 2 (R/lr™ > (2 - n) 

where c\, C2 and tt-hs are positive dimensionless constants. As in the St model, /(0) = 0, thus 
there is no formal cosmological constant. However, in high-curvature environments, namely 
when m 2 /R — > 0, an effective A-like term is present, since 

/ 2 \ n HS 

lim f(R)~-^m 2 + ^m 2 ( — ] . (2.12) 
m yR->o jy ' c 2 c 2 2 \R J V ' 

Finally, in the HS model the expansion history H in the early Universe is again the same as 
in the ACDM model. In the present case, the effective matter fraction is given by = 
6c 2 /(ci + 6c 2 ). 
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2.1 Background Evolution 

In a spatially flat Friedmann-Lemaitre-Robertson- Walker (FLRW) Universe with length ele- 
ment 

dx^dx u = -dt 2 + a 2 (t)dx 2 (2.13) 

filled with a perfect fluid with energy density p and pressure p, the modified Einstein's equa- 
tions (2.3) give a modified Friedmann's equation for the scale factor a = 1/(1 + z). By writing 
Eq. (2.3) in terms of the Hubble parameter H = dlna/dt, the Friedmann equation becomes 

H* + &H + tL*L* = *£ Pi (2.14) 
J,R VJ,R oj,R 

where a dot denotes a derivative with respect to the cosmic time t. We remind the reader 
that in a FLRW Universe, the Ricci tensor can be written as 

R = 6 ( H + 2H 



^H + 2H 2 y (2.15) 



Cardone et al. [55] solved the modified Friedmann equation (2.14) for a large set of 
randomly selected parameters {O m , X/c±, R*/c2, wst/^Hs} ; an d fit the expansion history 
of the Universe with a number of updated and accurate data. Specifically, they fitted the 
background expansion of the cosmos using the Hubble diagram of Snela [56] and Gamma 
Ray Bursts [55]. Moreover, they use H(z) data from passively evolving red galaxies [57], 
Baryon Acoustic Oscillations extracted from the seventh data release of the Sloan Digital Sky 
Survey, and distance priors from the recent WMAP7 data [58]. They found that the Hubble 
parameter can be excellently approximated by 

H(z) _ f Z(z)E CPL (z) + [1 - £(z)] E A (z) z < z A 

Ho ~ \ vfo$f>(i + *) s + i-n$f> z > z\ ( } 

where 

3 

£(z) = Y,ei{z-z A y (2.17) 

i=l 

is an interpolating function, with ej and z\ fitting parameters, and 

^CPL = n m (i + zf + (i - n m ) (i + zf^+^a) (218) 

is the dimensionless Hubble parameter in the phenomenological DE model of Chevalier, Po- 
larski and Linder [59, 60]. Eq. (2.18) reduces to E\ when the DE component is a cosmological 
constant A, i.e. (wq, w a ) = (—1, 0). The interpolating formula (2.16) of [55] means that the 
expansion rate H(z) for the two classes of f(R) models we are considering may be obtained 
by interpolating the CPL and the ACDM models back in time up to z\, whilst it becomes 
exactly the same as in ACDM at earlier times z > z\. In this ACDM- like era, the matter 
density of the model is an effective value given by 
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2.2 Evolution of Cosmological Perturbations 

In the linear theory of cosmological perturbations and in the Newtonian gauge, the metric 
(2.13) takes the form 

dx^dx u = - (1 + 2$) dt 2 + (1 + 2^) dx 2 , (2.20) 

with 

9nu = 9^ + ^9n^, (2-21) 

where ~g„ v is the background metric and 5g^ u the perturbation. The two scalar potentials 
<I> and \f/ are the metric perturbations in the Newtonian gauge. In this work we assume no 
anisotropic stress, thus the stress-energy tensor becomes 

Tf+v = T^ v + 5T^ U (2.22) 

with STq = —p5, 5T® = —pv yi and STj = 0, in the matter dominated epoch. Here, S = Sp/p 
is the density contrast, v is the velocity of the scalar perturbations and a comma denotes a 
derivative with respect to the spatial coordinates. 

Tsujikawa [61] showed that the evolution equation for the density contrast in f(R) gravity 
models is 

6 + 2HS - 4TT&(k, a)p5 ~ 0, (2.23) 
where there is an effective gravitational constant 

C 1 + 4^4^ 

g(M) = - - (2-24) 



f p 1 _j_ ok 2 J,RR 



The Poisson equation, which relates the potential $ to the distribution of matter over densities, 
in Fourier space reads 

$ fc (a) = -4tt^ (k,a)^pS k (a). (2.25) 

Finally, the other metric potential, is related to the above quantities by the parameter 
■q = — (<£> + which characterises the strength of the effective anisotropic stress. In the 

present case, it is given by 

2 k 2 f,RR 

V(k,a) = -^M-. (2.26) 



l+2 klLRR- 

a z f,R 



We remind that in standard GR, W(k, a) = G and rj = 0. 

The power spectrum of the matter fluctuations P s (k, z) is obtained by taking the Fourier 
transform of the two-point correlation function of the density contrast solution of Eq. (2.23). 
This means 

(6 k (z)5 k ,*(z)) = (2vr) 3 5 D (k - k') P 5 (k, z), (2.27) 
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where 5o is the Dirac 5 and k = |k|. Thus, the matter power spectrum can be written as 



P S (k,z) = 2tt 2 5 h 2 



T 2 (k) 



1 2 



5 k (z = 0) 



(2.28) 



with n s the tilt of the primordial power spectrum, 5h its normalisation and T{k) the matter 
transfer function, which describes the evolution of perturbations through the epochs of horizon 
crossing and radiation-matter transition. 

In the non-linear regime of the growth of matter overdensities, Eq. (2.23) does not hold 
any more. In the ACDM model, to obtain the nondinear matter power spectrum we need 
to resort to numerical solutions. A short cut is to use semi-analytical calculations [62] or 
fitting formulae extrapolated from numerical simulations (e.g. the halofit approach of [63]). 
For modified gravity to agree with Solar system observations, the non-linear matter power 
spectrum has to approach the standard ACDM solution on small scales. This means that 
the non-linear power spectrum has to be an interpolation of two power spectra. The former 
is the modified gravity non-linear power spectrum P^ /lG (k,z), which is obtained without the 
non-linear interactions that are responsible for the recovery of GR. This is equivalent to 
assume that gravity is modified down to small scales in the same way as in the linear regime. 
The latter term, P GR (k, z), is the non-linear power spectrum obtained in the DE model that 
follows the same expansion history of the Universe as the modified gravity model, yet obeying 
to GR. In other words, this is the non-linear power spectrum which will have a ACDM model 
with an expansions history H (z) equivalent to that of the modified gravity theory. 

Here, for the interpolation of these two spectra, we use the fitting formula proposed in 
Ref. [64] 

P 5 {h z) = P£ lG (k,z)+c Ql (z)tf(k,z)P^ R (k,z) ^ ^ 



l + c nl (z)^(k,z) 

where T, 2 (k,z) = [k 3 Pf in (k, z)/2ir 2 ] ai picks out non-linear scales, since Pf in (k,z) is the f(R) 
linear power spectrum, and c n \{z) = A (1 + z) a2 determines the scale at which the power 
spectrum approaches the GR result as a function of redshift. Their functional forms have 
been obtained by perturbation theory [65] and confirmed by A^-body simulations [66]. Here, 
we use a\ = 1/3, 02 = 1.05 and A = 0.08. 



3. Cosmic Shear Tomography 

In f(R) models, for the two potentials |$| ^ \^\ holds, even though there is no formal 
anisotropic stress in the energy momentum tensor T^ u . Thus, the relation between the distri- 
bution of matter overdensities in the Universe and the two metric perturbations, the potentials 
$ and ^, is not trivial. In GR, in the matter-dominated era, when there is no anisotropic 
stress, <3? = — \& and therefore we can simply use the Newtonian potential <£, thanks to the 
canonical Poisson equation, to compute cosmic convergence and shear. However, in general, 
the weak lensing effect is due to the combination of both the Newtonian and the metric po- 
tential. We will refer to this combination as the "deflecting potential," and we will denote it 
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with 2 

$ _ i|r 

T = —■ (3-1) 

Weak gravitational lensing is responsible for the shearing and magnification of the images 
of high-redshift sources due to the presence of intervening matter. The distortion are due 
to fluctuations in the gravitational potential, and are directly related to the distribution of 
matter and to the geometry and dynamics of the Universe. In particular, the distortions 
occurring to a background image can be decomposed into a convergence k and a (complex) 
shear 7 = 71 + ^72 , which are the entries of the distortion matrix 

f^( K + 71 72 V (3.2) 
\ 72 K-Ti/ 

The distortion matrix is directly related to background and perturbed cosmological quantities, 
since 

%i= [ X d X 'x'W(x')r^(n, X ') (3.3) 
Jo 

[44, 67], and commas denote derivatives with respect to directions perpendicular to the line 
of sight. Here, dx = dz/H(z) is the differential radial comoving distance, 

W( X ) = -2 X [°° d X ' ^W) (3.4) 

is the weight function of weak lensing, and n [x{z)\ represents the redshift distribution of the 
sources, such that fdxn(x) = 1. 

In the flat-sky approximation, we expand the shear 7(A) in its Fourier modes 

7 (n) = /(07We^. (3.5) 
The power spectrum is defined as the Fourier transform of the 2D correlation function 

(7(^)7*^)) = (^) 2 5 D (£ - £')C^£). (3.6) 

Thus, we have [42] 

where P r (k,z) is the power spectrum of the deflecting potential, and we have introduced 
Limber's approximation, where the only Fourier modes that contribute to the integral are 
those with £ = kx- 



2 In the literature, the combinations of the two potentials are often indicated by $± = — ($ ± ^) /2. 
However, different authors use to refer to the metric potentials in the length element differently. For the sake 
of simplicity, we have chosen this notation, because we are only interested in the combined effect responsible 
for the weak lensing signal. 
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In the case where one has distance information for individual sources, we can use this 
information for statistical studies. A natural course of action is to divide the survey into 
slices at different distances, and perform a study of the shear pattern on each slice [46]. In 
order to use the information effectively, it is necessary to look at cross-correlations of the 
shear fields in the slices, as well as correlations within each slice. This procedure is usually 
referred to as tomography. To better constrain these f(R) models with cosmic shear, we 
perform cosmic shear tomography. To do so, we separate the redshift distribution of sources 
n[x(z)] into redshift bins, each roughly containing the same number of sources. By doing so, 
the 2-dimensional shear power spectrum in the flat-sky approximation (3.7) reads 

^4/ d *™uv(|4 (38) 

where Wj(x) is the weight function (3.4) related to the ith bin. 



4. Fisher Matrix Analysis 

Cosmological parameters influence the shear in a number of ways, for instance the matter 
power spectrum P s (k,z) depends on Q m , h and the linear amplitude as- The linear power 
spectrum depends on the growth rate, which is also sensitive to the parameter of the A-like 
equation of state wa = pa/pa- It also affects the x( z ) relation and hence the angular diameter 
distance sinx [x( z )]- These parameters {t?a} may be estimated from the data using likelihood 
methods. Assuming uniform priors for the parameters, the maximum a posteriori probability 
for the parameters is given by the maximum likelihood solution. We use a Gaussian likelihood 

21nL = -Tr [in 6 + e -1 D] , (4.1) 

where C = ((d — d th )(d — d th ) T ) is the covariance matrix and D = (d — d th )(d — d th ) T is the 
data matrix, with d the data vector and d th the theoretical mean vector. 

The expected errors on the parameters can be estimated with the Fisher information 
matrix [68, 69, 70]. This has the advantage that different observational strategies can be 
analysed and this can be very valuable for experimental design. The Fisher matrix gives 
the best errors to expect, and should be accurate if the likelihood surface near the peak is 
adequately approximated by a multivariate Gaussian. 

The Fisher matrix is the expectation value of the second derivative of In L with respect 
to the parameters {# Q }, i.e. 

i 

and the marginal error on parameter # Q is 2 ■ If the means of the data are fixed, 

the Fisher matrix can be calculated from the covariance matrix and its derivatives [70] by 

? Q/3 = ^Tr [e-^aG-^p] . (4.3) 
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For a square patch of the sky, the Fourier transform leads to uncorrelated modes, provided 
the modes are separated by 27r/0 ra d where 6 ra( j is the side of the square in radians, and the 
Fisher matrix is simply the sum of the Fisher matrices of each £ mode, 

2£ + l 



(4.4) 



t 



where / s k y is the fraction of the sky covered by the survey under analysis, and 6 is the 
covariance matrix for a given £ mode. 

4.1 Bayesian Model Selection 

In this paper we compute parameter forecasts from cosmic shear tomography for viable f(R) 
models. It is worth noticing that we are dealing with a model alternative to standard ACDM. 
Hence, besides determining the best-fit value - and the errors - on a set of parameters within 
a model, we can also ask if this particular alternative model is preferable to the standard. 
Model selection is in a sense a higher-level question than parameter estimation. While in 
estimating parameters one assumes a theoretical model within which one interprets the data, 
in model selection one wants to know which theoretical framework is preferred given the 
data. Clearly, if our alternative model has more parameters than the standard one, chi- 
square analysis will not be of any use, because it will always reduce if we add more degrees of 
freedom. Bayesian analysis provides a useful "Occam's razor," which involves computation 
of the Bayesian evidence and of the Bayes factor B. 

The Bayesian evidence of a model M is defined as the marginalisation over the parameters 

p(d|M) = Jd m $p(d\d,M)p(<d\M), (4.5) 

where # is the parameter vector, and p(d|#,M) is the marginal likelihood of the parameters 
■d of the model M given the data d. Let us consider two competing models Mi and M2, the 
former nested in the latter. This means that Mi is simpler, because the set of its parameters 
is contained in the M2 parameter set {$ Q2 }, with a.\ = l,...,ni and Q2 = 1, . . . , ri2, 
ri2 > n\. In such a situation, one can compute the Bayes factor B, which is the ratio of 
the two corresponding posterior evidence probabilities p(Mi|d) and p(M2|d). The posterior 
probability for each model Mj is given by Bayes' theorem 

KM.|d) = "< d y . (4.6, 

If we choose noncommittal, constant priors p(M\) = p(M.2) = 1/2, the ratio of the 
posterior evidence probabilities reduces to the ratio of the evidence. Heavens et al. [71] 
showed that, in the Laplace approximation, where the expected likelihoods are given by 
multivariate Gaussians, and if one considers (B) as the ratio of the expected values, rather 
than the expectation value of the ratio, one eventually gets 

(B) = (2vr)-i f[ AW"*"*". (4-7) 

V 1 q=l 
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Here, 3~j is the Fisher matrix relative to the ith model, p = n<i — n\ is the number of extra 
parameters, and <5i9 is the vector of the parameter shifts. Indeed, if the correct underlying 
model is M 2 , the maximum of the expected likelihood will not, in general, be at the correct 
parameter values of Mi (see Fig. 1 of [71]). The n\ parameters of Mi shift their values to 
compensate the fact that i9 ai +i, . . . , flai+p are kept fixed at an incorrect fiducial value $ai+i = 
■ ■ ■ = $ ai +p = 0. These shifts can be computed under the assumption of a multivariate 
Gaussian distribution [72], and read 

s& = -9 r r 1 g 2 ^, (4.8) 

with S2 a subset of the M2 Fisher matrix and 8ij) the shifts of the p extra parameters ip. 

It is usual to consider the logarithm of the Bayes factor, for which the so-called "Jeffreys' 
scale" gives empirically calibrated levels of significance for the strength of evidence [73]. A 
more recent version of Jeffreys' scale sets 1 < | lnB\ < 2.5 as "substantial" evidence in favour 
of a model, 2.5 < \ lnB\ < 5 as "strong," and |lnB| > 5 as "decisive." These descriptions 
seem too aggressive, for | kiB\ = 1 corresponds to a posterior probability for the less-favoured 
model which is 0.37 times the favoured model [74]. Other authors have introduced different 
terminology (e.g. [75]). 

5. Results 

We compute our results for a flat Universe with cosmological parameters h = 0.7, Sl m = 0.28, 

n b = 2.22 • w- 2 h- 2 , n A = i- n m , i ogl0 a/ i ogl0 Cl = 2.38/4.98, io gl0 «/ i ogl0 c 2 = -2.6/3.79 

and nst/nus = 1.79/1.64, where k = R*/Ro, as obtained with the Monte Carlo Markov 
Chains of Ref. [55]. For the matter power spectrum, we use the transfer function proposed 
by Eisenstein & Hu (1998) [76], rms mass fluctuations as = 0.8 on a scale of Sh~ l Mpc, and 
spectral index n s = 0.96 [58]. 

In f(R) models, gravity is stronger than in GR [61]. In particular, the effective gravi- 
tational constant &(k,a), which appears in the source term driving the evolution of matter 
density perturbations (Eq. 2.23), can change significantly compared to the Newtonian gravi- 
tational constant G. It has been shown that in the so-called "scalar-tensor" regime & ~ 4G/3 
[61]. In Fig. 1 we show the present-day value ^(/c,a = ao) of the effective gravitational con- 
stant normalised to the Newtonian G. At large scales, f(R) gravity behaves like GR, whilst 
the larger is the value of k, the greater is the difference. In more detail, the HS model presents 
a gravitational coupling larger than the St model, because it reaches the 4G/3 ~ 1.33G value 
at smaller scales. 

5.1 The Cosmic Shear Signal 

We compute our results for a 20, 000 deg 2 next-generation cosmic shear experiment such as 
Euclid. 3 We use the Euclid baseline given in the Yellow Book [77]. The source distribution 

3 http://sci.esa.int/science-e/www/area/index.cfm?fareaid=102 
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Figure 1: The present-day effective gravitationai constant 
scale fc, normalised to the Newtonian constant G. 



&(k 7 a = a ) versus the physical 



over redshifts has the form [78] 



n{z) oc z 2 e 



-(*)' 



(5.1) 



where zq = z m /1.4, and z m = 0.89 is the median redshift of the survey. The number density of 
the sources, with estimated redshift and shape, is 35 per square arcminute. We also compute 
the expected errors according to [42, 44] 



AC^(£) = 



(2£+l)f sky 



CT{1) + 



(Tint 2 ) 1 



n 



(5.2) 



where / s k y = 6^^/129, 600 is the fraction of the sky covered by a survey of area ©deg an< ^ 
(Tint 2 ) ' 5 — 0.4 is the galaxy-intrinsic shear rms in one component. 



The peculiarities of the two f(R) models depicted in Fig. 1 propagate into the cosmic 
shear power spectrum C(£), which is shown in Fig. 2. The curves refer to the standard 
ACDM cosmology (solid), the St (dashed) and the HS (dot-dashed) models, respectively. 
At large angular scales the three signals are indistinguishable, while for large values of £, 
differences start to be important. This is a consequence of Limber's approximation, which 
sets k = £/x and therefore shows at large £'s the features %(/c) presents at large fc's. Thus, we 
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observe in the cosmic shear power spectrum the same behaviour of the effective gravitational 
constant a). Specifically, the HS signal differs from ACDM by more than la at i > 300. 
On the contrary, the St model produces a cosmic shear power spectrum in agreement with 
ACDM up to t ~ 5000, and the errorbars show that it is almost completely degenerate with 
ACDM. 
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Figure 2: Shear power spectrum £(£ + 1)C~> (£) / {2n) in the ACDM (solid), St (dashed) and HS 
(dot-dashed) models; la errorbars on the ACDM signal, computed according to Eq. (5.2), are also 
shown. 

To better constrain the model parameters, and with the aim of lifting the degeneracy 
between ACDM and the St model, we calculate the tomographic shear matrix Cj-{t). To do 
so, we separate the Euclid distribution of sources (5.1) into ten redshift bins. Fig. 3 presents 
the diagonal elements Cj i {£) of the tomographic shear matrix as a function of the angular scale 
£ for standard ACDM cosmology (solid), the St (dashed) and the HS (dot-dashed) models, 
respectively. We also show la errorbars on the ACDM signal. Though tomography does 
enhance the quality of the signal when one looks at the high-redshift autocorrelations, we 
find that the degeneracy between ACDM and the St model is not removed. The off-diagonal 
elements of the tomographic shear matrix show a similar behaviour. In principle, they should 
be more useful, since the Poissonian noise term in Eq. (5.2) holds for correlations between 
the same bin only. However, the St signal is still too close to what predicted by ACDM, at 
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Figure 3: Diagonal elements C^(£) of the tomographic shear matrix for the ACDM (solid), St (dashed) 
and HS (dot-dashed) models. Curves refer to shear autocorrelations, from the l st -l st to the 10 th -10 th 
bin pairs with increasing redshift; la errorbars on the ACDM signal, computed according to Eq. (5.2), 
are also shown. 

least for the range of angular scales probed here. 
5.2 Cosmological Parameter Forecasts 

Once we have introduced the method and the survey design formalism, we now show the 
cosmological parameter forecasts for such a survey. By using the Fisher matrix analysis 
outlined in Ref. [72] , we calculate the predicted Fisher matrices and the parameter constraints 
for a 20,000deg 2 Euclid-like survey. In all Fisher matrix calculations we use a six-parameter 
cosmological set {Q m = ft D M + ^6, h, log 10 A/log 10 c±, log 10 k/ log 10 c 2 , n st /n H s}- Note 
that we do not use Qa as a free parameter. 

It is worth noting that, in the high-wavenumber regime, the fitting formulae of [63] may 
be unreliable, or baryonic effects might alter the power spectrum (k > 10/iMpc" 1 [79, 80]). 
Hence, we do not analyse modes with k > 1.5/iMpc _1 . However, in Limber's approximation 
I oc A; -1 , the choice of the maximum angular wavenumber ^ max entering Eq. (4.4) is crucial. 
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Table 1: Forecasts on the St model parameters. 



The largest proper £ max is a matter of debate (e.g. [81, 82]). Here, we follow three approaches, 
one more conservative with £ max = 1000, one optimistic with £ max = 5000 - which was already 
used in testing alternative cosmological models with cosmic shear [83] -, and one where £ max 
is determined in a bin-dependent way. Actually, fixing the distances of the tomographic 
binning removes some flexibility in the probed physical wavenumbers; so there is a risk that 
some useful modes are excluded (thus increasing the statistical errors), and/or that, for the 
nearby shells, the sampled physical wavenumber range extends to too high a value of k, 
where theoretical uncertainties become a potential source of systematic error. To avoid these 
potential problems, we define £ max = x( z bin)/k max [84]: we thus increase the largest angular 
scale for distant shells and reduce it for nearby shells. 

Figs. 4-5 show the la and 2a contours obtained by calculating the Fisher matrices of 
the tomographic shear signal (3.8), for the St and HS models, respectively. The green curves 
refer to £ max = 1000, the blue ones to £ max (z\,i n ) and the red ones to £ max = 5000. 

The use of tomography enhances the accuracy significantly, compared to the cosmic 
shear power spectrum. This analysis thus yields tighter constraints on the model parameters. 
Specifically, we show our constraints in Tables 1-2, for the St and HS models, respectively, 
where ^ is the mean fiducial value of the parameter and a$ is the predicted standard de- 
viation. However, these small constraints may probably under-predict the true parameter 
errors. This is due to the fact that gravitational lensing does mix angular modes to some 
degree because of its intrinsic non-local structure. Non-linear mode coupling affects Fisher 
matrices for cosmological observations. Indeed, these couplings tend to correlate small-scale 
power, moving information from lower to higher-order moments of the shear field. Therefore, 
Gaussian approximations may produce over-optimistic forecasts [85]. 

5.3 Model Selection 

In Section 4.1 we showed how the Bayes factor can be used to determine which model is 
favoured by the data. We outlined how the issue of Bayesian model selection can be performed 
in the case of nested models. Actually, only the HS model is formally nested in ACDM. Indeed, 
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Figure 4: Expected marginal errors on the St model cosmological parameters from a 20, 000 deg 2 
Euclid-like survey. Ellipses show the la and 2a errors for two parameters (68.3% and 95.8% confidence 
regions, respectively), marginalised over all the other parameters. The green, red and blue curves refer 
to the conservative, bin-dependent and optimistic £ max . 

if one sets nns = 0, the two extra parameters, c\ and C2, depend on each other. Specifically, 

log loCl = 6(l + 10 10 ^)^ (53) 

must hold. 

By using the Fisher matrix formalism for a Euclid-like survey, we compute the Bayes 
factor In B as a function of nns for the HS model over the standard ACDM cosmology. The 
result is presented in Fig. 6. The green, red and blue curves refer to £ max = 1000, 5000 and 
Anax^bin)] respectively. A solid line denotes a positive value of In B, whereas a dashed line 
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Figure 5: Same as Fig. 4 for the HS model parameters. 
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Table 2: Same as Table 1 for the HS model. 
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means \n.B < 0. Therefore, if Euclid detects any nonzero nns's, the Bayes factor will favour 
the HS model only on the right of the cusps. Contrarily, Occam's razor will prefer ACDM 
anyway. The three dot-dashed horizontal lines show the "strength" of this favour /disfavour, 
being them InB = 1, 2.5 and 5, as in the current version of Jeffreys' scale. We want to 
emphasise that the HS model value uhs = 1-64 lies in the region of "decisive" evidence in 
favour of the more complex model. Thus, Euclid will favour the HS model unquestionably, if 
nns = 1-64 will be measured. 




0.01 0.1 1 

^HS 

Figure 6: The Bayes factor InB for the HS model (nested case) over standard ACDM as a function 
of the extra parameter nus- The green, red and blue curves refer to the conservative, bin-dependent 
and optimistic £ max , respectively. The horizontal lines denote the Jeffreys' scale levels of significance. 

6. Conclusions 

In this paper, we study whether the cosmic shear signal measured by the ESA Cosmic Vision 
Euclid satellite [86, 77, 87] can constrain two viable f(R) models. These models, proposed 
by Starobinsky [40] (St) and Hu & Sawicki [41] (HS) to be able to pass Solar System tests, 
account for the present-day accelerated expansion of the Universe without any cosmological 
constant term. At the same time, the functional form of their f(R) correctly reproduces the 
Hubble rate of an early-time matter-dominated cosmos which undergoes a late-time phase of 
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accelerated expansion. Recently, Thomas et al. [88] performed a similar analysis, but they 
focused on constraining the scalaron mass and its parameterisations. 

Cardone et al. [55] derived the St and HS model parameters against several cosmological 
datasets by fitting the expansion rate history H(z). Specifically, they used the Hubble diagram 
of Snela [56] and Gamma Ray Bursts [55]. Moreover, they use H(z) data from passively 
evolving red galaxies [57], Baryon Acoustic Oscillations extracted from the seventh data 
release of the Sloan Digital Sky Survey, and distance priors from the recent WMAP7 data 
[58]. Therefore, we use Cardone et al. [55] values for the model parameters to compute our 
observables. 

We calculate the cosmic shear power spectra and tomographic shear matrices for the two 
models and compared them with ACDM. Beynon et al. [89] showed that there is substantial 
additional discriminatory power between f(R) models by the inclusion of the non-linear power 
regime. They have also shown that using only the halofit formula [63] without any attempt 
to obtain the GR non-linear power spectrum on small scales leads to an overestimation in the 
ability of future surveys to differentiate between different growth histories. Similarly, Casarini 
et al. [90] demonstrated that halofit expressions, well tested for the ACDM model, implies 
substantial discrepancies with respect to results directly obtained from A-body simulations, 
when the effective wbe(z) 7^—1- To avoid all these problems, we use the fitting formulae 
proposed in Ref. [64], which are an interpolation of the modified-gravity non-linear matter 
power spectrum and the small-scales GR prediction. These formulas have been confirmed 
both theoretically [65] and by A-body simulations [66], although for a slightly different f(R) 
model. 

In f{R) models, gravity is stronger than in GR [61], and in particular the Newtonian 
constant of gravitation is replaced by a time- and scale-dependent effective &(k,a). Further- 
more, when gravity is not described by GR, the two metric potentials $ and ^ are no longer 
equal in modulus. Therefore, one has to take the deflecting potential Y = — (3> — 1 I r )/2 into 
account, when computing lensing observables. As expected, we find that the cosmic shear 
signal tracks the differences between the effective gravitational constant in the St, HS and 
ACDM models. Specifically, the HS shear signal is different from ACDM at angular scales 
£ > 300, whilst for the St model the agreement holds for all i up to £ ~ 5000. We also find 
that the difference between the ACDM cosmic shear signal and what is expected for the St 
model are too weak to allow one to discriminate between them. 

Then, we exploit the power of cosmic shear tomography to constrain the model param- 
eters. We separate the Euclid distribution of sources into ten redshift bins, each containing 
roughly the same number of sources. In the present study, we neglect the contribution of 
intrinsic alignments of the source shapes. We find that Euclid is able to tightly constrain 
the parameters of these f(R) models, as shown in Tables 1-2. Since the choice of the largest 
angular scale entering in the computation of the Fisher matrix is of primary importance, we 
perform the Fisher analysis for more f max 's. In particular, we use a conservative value of 
1000, an optimistic value of 5000 and a bin-dependent setting, which increases the maximum 
angular wavenumber for distant shells and reduces it for nearby shells. 
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Finally, we compute the Bayesian expected evidence (e.g. [75, 71]) for the HS model over 
the ACDM model as a function of the extra parameter nns- This can be done because the 
ACDM model is formally nested in the HS model, and the latter is equivalent to the former 
when h-hs = 0. In this case, the two other extra parameters, c\ and C2, are linked together 
by the request of reproducing a A-like term in the gravitational Lagrangian. The expected 
evidence clearly shows that the Euclid survey data will unquestionably favour the HS model 
if any value nns > 0.02 is measured - according to either the optimistic or the bin-depending 
analysis. 
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